*********************************************************************************
* Table 7. Beliefs – IPWRA (analytic SEs, main text) + OLS (appendix)
*********************************************************************************

	****************************************************
	* 	SET 1
	
	*	yvars = fob2   sob2p Dmis_sob2p  fob3  sob3p Dmis_sob3p
	
	*	sample3_ext & sample2 & sample4_any
	****************************************************


* Main sample indicator
local samples  sample4_other /*sample3_ext sample2 sample4_any  sample4_both  sample4_other*/ /* sample4_alone*/

* Baseline belief block (used in diagnostics)
local bel2  fob2 sob2 sob2m sob2p

* "Strata" covariates used in absorb() and IPW models
local strat inactiveF employedF fob2M

* Outcomes: comment/uncomment to include/exclude
	local yvars fob2   sob2p Dmis_sob2p  fob3  sob3p Dmis_sob3p
	*local yvars2 fob2 DDisagree_fob2p D  DExpDisagree_sob2p_w3  fob3 DDisagree_fob3p DExpDisagree_sob3p_w3

* Initial control variables (diagnostics)
local xvars female c1 c2_2_rp c3 c5 f6 age_dif  level_1 level_2 edu_WmoreH  inactiveM unemployedM estrato_12 estrato_3 tec_2 


* Belief blocks
local bel1 fob1 sob1 sob1sdb sob1m 
local bel2 fob2 sob2 sob2m sob2p
local bel3 fob3 sob3 sob3sdb sob3m sob3p
local bel4 fob4 fob6
local bel5 fob5 sob5 sob5m sob5p
local bel6 fob7 sob7 sob7m sob7p
local bel7 fob8 sob8 sob8m sob8p
*local bel8 e501 e502 e503 e501b e504 e505 e506 e504b
local bel9 fob9

* All belief variables used in selection model
local bels `bel1' `bel2' `bel3' `bel4' `bel6' `bel7'  `bel9'

		foreach sample in `samples' {

********************************************************
* Diagnostics for attrition (who is in sample3_ext)
********************************************************

reghdfe `sample' D `xvars' `bel2', absorb(`strat') vce(cluster id)
	/* sign : D, female, age, level_1, inactivM, tec_2 */
reg sample2 sample1_ext `xvars' `bel2' `strat' if D==0,   vce(cluster id)
	/* sign: being in whatsapp; (W more edu H) ; sob2p; tec_2 */	
reg sample2 sample1_ext `xvars' `bel2' `strat' if D==1,   vce(cluster id)
	/* sign: being in whatsapp;female; age; age_dif;  tec_2; sob2m: inactiveF */

***********	
reg sample2 sample1_ext##D `xvars' `bel2' `strat' ,   vce(cluster id)
	/* sign: being in whatsapp;female; age;  tec_2 */
***********		
	
reg sample1_ext D `xvars' `bel2' `strat' ,   vce(cluster id)
	/* sign: female c5 age, age dif, (level_1), inactiveM, tec_2 */

*********************************************************	
* Diagnostics for selective exposure (who takes D)
*********************************************************

reghdfe D `xvars' `bel2' if `sample' == 1, absorb(`strat') vce(cluster id)

reg D `xvars' `bel2' `strat' if sample3 == 1,   vce(cluster id)
	/* sign: sob2p */
	
****************************************************
* According to diagnostics, our working controls:
* c5 level_1 tec_2 fob2 sob2p
* Working control set for main outcome models
local vars_sample3 `xvars' `bel2'

********************************************************************************
* IPW weights for selection (attrition weights for being in `sample')
********************************************************************************
ta sample3 D ,m
ta sample3 D if sample3_ext ==1,m

cap drop ps_attrition 
cap drop w_*
cap drop wa_*
cap drop ws_*

	probit `sample' D `xvars' `bels' `strat', vce(cluster id)
	*probit `sample' D `xvars' `bels' `strat', vce(cluster id)

predict ps_attrition, pr 

gen w_`sample' = 1 / ps_attrition
quietly summ w_`sample'
local p_uncond = r(mean)
gen ws_`sample' = `p_uncond' / ps_attrition   // normalized (if needed)
gen wa_`sample' = 1                           // unweighted

* Default selection weight (you can switch to ws_ if you prefer normalized)
local weight w_`sample'
summ `weight' 
bys D : su `weight' if sample3_ext ==1
bys D : su `weight' ws_`sample' ps_attrition if sample3_ext ==1
local reps $reps

********************************************************************************
* STORAGE MACROS FOR THE THREE TABLES
********************************************************************************
eststo clear
local j = 1

local estlist_ols_base ""   // 7A: Baseline OLS (no weights) – one model per outcome
local mtitles_ols_base ""

local estlist_ols_ipw ""    // 7B: OLS + selection IPW – one model per outcome
local mtitles_ols_ipw ""

local estlist_ipwra ""      // 7C: IPWRA ATT/PO + OLS-IPW D – one model per outcome
local mtitles_ipwra ""


********************************************************************************
* LOOP OVER OUTCOMES
********************************************************************************

foreach y of local yvars{

    /********************************************************************
     * 1. BASELINE OLS (no selection weights): pooled model only
     *    Male/Female enter as stats rows
     ********************************************************************/

    *** (1.1) Pooled OLS – All (unweighted)
    reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1, absorb(`strat') vce(cluster id)
    scalar N_all = e(N)
    eststo OLS_base_`j'

    *** Means of Y by subgroup (unweighted)
    quietly summarize `y'_w3 if `sample' == 1
    scalar ymean_all  = r(mean)

    quietly summarize `y'_w3 if `sample' == 1 & female == 0
    scalar ymean_male = r(mean)

    quietly summarize `y'_w3 if `sample' == 1 & female == 1
    scalar ymean_fem  = r(mean)

    *** (1.2) OLS – Male only (unweighted)
    reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1 & female == 0, ///
        absorb(`strat') vce(cluster id)
    scalar N_male     = e(N)
    scalar b_D_male   = _b[D]
    scalar se_D_male  = _se[D]
    scalar t_D_male   = b_D_male / se_D_male
    scalar p_D_male   = 2*(1 - normal(abs(t_D_male)))

    *** (1.3) OLS – Female only (unweighted)
    reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1 & female == 1, ///
        absorb(`strat') vce(cluster id)
    scalar N_fem      = e(N)
    scalar b_D_fem    = _b[D]
    scalar se_D_fem   = _se[D]
    scalar t_D_fem    = b_D_fem / se_D_fem
    scalar p_D_fem    = 2*(1 - normal(abs(t_D_fem)))

    *** RITEST p-values for D (All, Male, Female), with capture to avoid r(111)
    * All
    capture noisily ritest D _b[D], ///
        cluster(id) strata(`strat') reps(`reps') : ///
        reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1, ///
        absorb(`strat') vce(cluster id)
    scalar p_ri_all = .
    if (_rc == 0) {
        matrix P_all = r(p)
        scalar p_ri_all = round(P_all[1,1], .001)
    }

    * Male
    capture noisily ritest D _b[D], ///
        cluster(id) strata(`strat') reps(`reps') : ///
        reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1 & female == 0, ///
        absorb(`strat') vce(cluster id)
    scalar p_ri_male = .
    if (_rc == 0) {
        matrix P_m = r(p)
        scalar p_ri_male = round(P_m[1,1], .001)
    }

    * Female
    capture noisily ritest D _b[D], ///
        cluster(id) strata(`strat') reps(`reps') : ///
        reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1 & female == 1, ///
        absorb(`strat') vce(cluster id)
    scalar p_ri_fem = .
    if (_rc == 0) {
        matrix P_f = r(p)
        scalar p_ri_fem = round(P_f[1,1], .001)
    }

    *** Attach all scalars to the pooled OLS model for 6A
    estimates restore OLS_base_`j'
    estadd scalar N_all      = N_all,      replace
    estadd scalar N_male     = N_male,     replace
    estadd scalar N_fem      = N_fem,      replace

    estadd scalar ymean_all  = ymean_all,  replace
    estadd scalar ymean_male = ymean_male, replace
    estadd scalar ymean_fem  = ymean_fem,  replace

    estadd scalar b_D_male   = b_D_male,   replace
    estadd scalar se_D_male  = se_D_male,  replace
    estadd scalar p_D_male   = p_D_male,   replace

    estadd scalar b_D_fem    = b_D_fem,    replace
    estadd scalar se_D_fem   = se_D_fem,   replace
    estadd scalar p_D_fem    = p_D_fem,    replace

    estadd scalar p_ri_all   = p_ri_all,   replace
    estadd scalar p_ri_male  = p_ri_male,  replace
    estadd scalar p_ri_fem   = p_ri_fem,   replace

    eststo OLS_base_`j'

    /********************************************************************
     * 2. OLS WITH SELECTION IPW WEIGHTS: pooled model only
     ********************************************************************/

    *** (2.1) Pooled weighted OLS – All
    reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1 [pweight = `weight'], ///
        absorb(`strat') vce(cluster id)
    scalar N_all_w = e(N)
    eststo OLS_ipw_`j'

    *** Weighted means of Y by subgroup
    quietly summarize `y'_w3 [aweight = `weight'] if `sample' == 1
    scalar ymean_all_w  = r(mean)

    quietly summarize `y'_w3 [aweight = `weight'] if `sample' == 1 & female == 0
    scalar ymean_male_w = r(mean)

    quietly summarize `y'_w3 [aweight = `weight'] if `sample' == 1 & female == 1
    scalar ymean_fem_w  = r(mean)

    *** (2.2) Weighted OLS – Male
    reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1 & female == 0 [pweight = `weight'], ///
        absorb(`strat') vce(cluster id)
    scalar N_male_w     = e(N)
    scalar b_D_male_w   = _b[D]
    scalar se_D_male_w  = _se[D]
    scalar t_D_male_w   = b_D_male_w / se_D_male_w
    scalar p_D_male_w   = 2*(1 - normal(abs(t_D_male_w)))

    *** (2.3) Weighted OLS – Female
    reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1 & female == 1 [pweight = `weight'], ///
        absorb(`strat') vce(cluster id)
    scalar N_fem_w      = e(N)
    scalar b_D_fem_w    = _b[D]
    scalar se_D_fem_w   = _se[D]
    scalar t_D_fem_w    = b_D_fem_w / se_D_fem_w
    scalar p_D_fem_w    = 2*(1 - normal(abs(t_D_fem_w)))

    *** Attach scalars to pooled weighted OLS model for 6B
    estimates restore OLS_ipw_`j'
	estadd scalar N_all      = N_all_w,      replace
    estadd scalar N_male     = N_male_w,     replace
    estadd scalar N_fem      = N_fem_w,      replace

    estadd scalar ymean_all  = ymean_all_w,  replace
    estadd scalar ymean_male = ymean_male_w, replace
    estadd scalar ymean_fem  = ymean_fem_w,  replace

    estadd scalar b_D_male   = b_D_male_w,   replace
    estadd scalar se_D_male  = se_D_male_w,  replace
    estadd scalar p_D_male   = p_D_male_w,   replace

    estadd scalar b_D_fem    = b_D_fem_w,    replace
    estadd scalar se_D_fem   = se_D_fem_w,   replace
    estadd scalar p_D_fem    = p_D_fem_w,    replace

    eststo OLS_ipw_`j'

    /********************************************************************
     * 3. IPWRA – ATT/PO for All, Male, Female (analytic SEs)
     *    Use selection weights as pweights; attach to pooled weighted model
     ********************************************************************/

    *** (3.1) IPWRA – All
    teffects ipwra (`y'_w3 `y' `strat') ///
                   (D `vars_sample3' `strat') ///
                   if `sample' == 1 [pweight = `weight'], ///
                   atet vce(cluster id)

	scalar N_all_ipw = e(N)   
    matrix bA = e(b)
    matrix VA = e(V)

    scalar ATET_all    = bA[1,1]
    scalar ATET_se_all = sqrt(VA[1,1])
    scalar ATET_pv_all = 2*(1 - normal(abs(ATET_all/ATET_se_all)))

    scalar PO_all      = bA[1,2]
    scalar PO_se_all   = sqrt(VA[2,2])
    scalar PO_pv_all   = 2*(1 - normal(abs(PO_all/PO_se_all)))

    *** (3.2) IPWRA – Male
    teffects ipwra (`y'_w3 `y' `strat') ///
                   (D `vars_sample3' `strat') ///
                   if `sample' == 1 & female == 0 [pweight = `weight'], ///
                   atet vce(cluster id)

    scalar N_male_ipw = e(N) 
	matrix bA = e(b)
    matrix VA = e(V)

    scalar ATET_male    = bA[1,1]
    scalar ATET_se_male = sqrt(VA[1,1])
    scalar ATET_pv_male = 2*(1 - normal(abs(ATET_male/ATET_se_male)))

    scalar PO_male      = bA[1,2]
    scalar PO_se_male   = sqrt(VA[2,2])
    scalar PO_pv_male   = 2*(1 - normal(abs(PO_male/PO_se_male)))

    *** (3.3) IPWRA – Female
    teffects ipwra (`y'_w3 `y' `strat') ///
                   (D `vars_sample3' `strat') ///
                   if `sample' == 1 & female == 1 [pweight = `weight'], ///
                   atet vce(cluster id)
	
	scalar N_fem_ipw = e(N)
    matrix bA = e(b)
    matrix VA = e(V)

    scalar ATET_fem    = bA[1,1]
    scalar ATET_se_fem = sqrt(VA[1,1])
    scalar ATET_pv_fem = 2*(1 - normal(abs(ATET_fem/ATET_se_fem)))

    scalar PO_fem      = bA[1,2]
    scalar PO_se_fem   = sqrt(VA[2,2])
    scalar PO_pv_fem   = 2*(1 - normal(abs(PO_fem/PO_se_fem)))

    *** Create IPWRA container: copy pooled weighted OLS and attach ATT/PO scalars
    estimates restore OLS_ipw_`j'
    eststo IPWRA_`j'

    estimates restore IPWRA_`j'
	
	estadd scalar N_all_ipw  = N_all_ipw,  replace
    estadd scalar N_male_ipw = N_male_ipw, replace
    estadd scalar N_fem_ipw  = N_fem_ipw,  replace
	
	
    estadd scalar ATET_all    = ATET_all,    replace
    estadd scalar ATET_se_all = ATET_se_all, replace
    estadd scalar ATET_pv_all = ATET_pv_all, replace

    estadd scalar PO_all      = PO_all,      replace
    estadd scalar PO_se_all   = PO_se_all,   replace
    estadd scalar PO_pv_all   = PO_pv_all,   replace

    estadd scalar ATET_male     = ATET_male,     replace
    estadd scalar ATET_se_male  = ATET_se_male,  replace
    estadd scalar ATET_pv_male  = ATET_pv_male,  replace

    estadd scalar PO_male       = PO_male,       replace
    estadd scalar PO_se_male    = PO_se_male,    replace
    estadd scalar PO_pv_male    = PO_pv_male,    replace

    estadd scalar ATET_fem      = ATET_fem,      replace
    estadd scalar ATET_se_fem   = ATET_se_fem,   replace
    estadd scalar ATET_pv_fem   = ATET_pv_fem,   replace

    estadd scalar PO_fem        = PO_fem,        replace
    estadd scalar PO_se_fem     = PO_se_fem,     replace
    estadd scalar PO_pv_fem     = PO_pv_fem,     replace

    eststo IPWRA_`j'

    ********************************************************************
    * Append models and column titles
    ********************************************************************
    local ylab : variable label `y'_w3
    if "`ylab'" == "" local ylab "`y'_w3"

    local estlist_ols_base "`estlist_ols_base' OLS_base_`j'"
    local mtitles_ols_base `"`mtitles_ols_base' "`ylab'""'

    local estlist_ols_ipw  "`estlist_ols_ipw'  OLS_ipw_`j'"
    local mtitles_ols_ipw  `"`mtitles_ols_ipw'  "`ylab'""'

    local estlist_ipwra    "`estlist_ipwra'    IPWRA_`j'"
    local mtitles_ipwra    `"`mtitles_ipwra'    "`ylab'""'

    local j = `j' + 1
}

/*
*=====================================================
* Table 7A – BASELINE OLS (no selection weights)
*=====================================================
#delimit ;
local opts_ols_base style(tex)  
  varlabels(_cons Constant D "Treatment" female "Female (=1)") 
  cells(b(star fmt(%9.3f)) se(par)) 
  stats(N_all N_male N_fem 
        ymean_all ymean_male ymean_fem
        b_D_male se_D_male p_D_male
        b_D_fem  se_D_fem  p_D_fem
        p_ri_all p_ri_male p_ri_fem,
        fmt(%9.3f)
        label("N (All)" 
              "N (Male)"
              "N (Female)"
              "Mean Dep. Var (All)"
              "Mean Dep. Var (Male)"
              "Mean Dep. Var (Female)"
              "OLS β(D) Male"    "OLS SE Male"    "OLS p-val Male"
              "OLS β(D) Female"  "OLS SE Female"  "OLS p-val Female"
              "RI p-val D (All)" "RI p-val D (Male)" "RI p-val D (Female)")) 
  label collabels(none) msign(--) nolz varwidth(1) modelwidth(14) 
  starlevels(* 0.10 ** 0.05 *** 0.01)
;
#delimit cr

estout `estlist_ols_base' using "$outdir/Table7_SpousalBeliefs_OLS_baseline_`sample'.tex", ///
    replace `opts_ols_base' mlabels(`mtitles_ols_base')

*=====================================================
* Table 7B – OLS WITH SELECTION IPW WEIGHTS
*=====================================================
#delimit ;
local opts_ols_ipw style(tex)  
  varlabels(_cons Constant D "Treatment" female "Female (=1)") 
  cells(b(star fmt(%9.3f)) se(par)) 
  stats(N_all N_male N_fem 
        ymean_all ymean_male ymean_fem
        b_D_male se_D_male p_D_male
        b_D_fem  se_D_fem  p_D_fem,
        fmt(%9.3f)
        label("N (All)" 
              "N (Male)"
              "N (Female)"
              "Mean Dep. Var (All, weighted)"
              "Mean Dep. Var (Male, weighted)"
              "Mean Dep. Var (Female, weighted)"
              "OLS-IPW β(D) Male"    "OLS-IPW SE Male"    "OLS-IPW p-val Male"
              "OLS-IPW β(D) Female"  "OLS-IPW SE Female"  "OLS-IPW p-val Female")) 
  label collabels(none) msign(--) nolz varwidth(1) modelwidth(14) 
  starlevels(* 0.10 ** 0.05 *** 0.01)
;
#delimit cr

estout `estlist_ols_ipw' using "$outdir/Table7_SpousalBeliefs_OLS_IPW_`sample'.tex", ///
    replace `opts_ols_ipw' mlabels(`mtitles_ols_ipw')
*/
*=====================================================
* Table 7C – IPWRA (analytic SEs), only D in coefficient panel
*=====================================================
#delimit ;
local opts_ipwra style(tex)  
  cells(b(star fmt(%9.3f)) se(par)) 
  stats(ATET_all ATET_se_all ATET_pv_all 
        PO_all   N_all_ipw
        ATET_male ATET_se_male ATET_pv_male
        PO_male N_male_ipw
		ATET_fem  ATET_se_fem  ATET_pv_fem
        PO_fem    N_fem_ipw,
        fmt(%9.3f)
        label("\textbf{ATT All}"      "\hspace{3mm}\textit{(Std. Err)}"   "P-value ATT All"
              "\hspace{3mm}\textit{Mean Dep.\ Var (Controls)}"       "\hspace{3mm} N All"
              "\textbf{ATT Male}"     "\hspace{3mm}\textit{(Std. Err)}"  "P-value ATT Male"
              "\hspace{3mm}\textit{Mean Dep.\ Var (Controls)}"      "\hspace{3mm} N Male"
              "\textbf{ATT Female}"   "\hspace{3mm}\textit{(Std. Err)}""P-value ATT Female"
              "\hspace{3mm}\textit{Mean Dep.\ Var (Controls)}"    "\hspace{3mm}N Female")) 
  label collabels(none) msign(--) nolz varwidth(1) modelwidth(14) 
  starlevels(* 0.10 ** 0.05 *** 0.01)
;
#delimit cr

estout `estlist_ipwra' using "$outdir/Table7_SpousalBeliefs_IPWRA_main_`sample'.tex", ///
    replace `opts_ipwra' mlabels(`mtitles_ipwra') ///
    keep(D) varlabels(D "Treatment") 

}

/*
*********************************************************************************************

*********************************************************************************
* Table 7. Beliefs – IPWRA (analytic SEs, main text) + OLS (appendix)
*********************************************************************************

	****************************************************
	* 	SET 2
		
	*	yvars = DDisagree_fob2p DExpDisagree_sob2p_w3 Dmis_sob2p    DDisagree_fob3p DExpDisagree_sob3p_w3 Dmis_sob3p

	
	****************************************************


* Main sample indicator
local samples  sample3_ext sample4_any /*sample2 sample4_other  sample4_both sample4_other sample4_alone*/

* Baseline belief block (used in diagnostics)
local bel2  fob2 sob2 sob2m sob2p

* "Strata" covariates used in absorb() and IPW models
local strat inactiveF employedF fob2M

* Outcomes: comment/uncomment to include/exclude
	*local yvars fob2   sob2p Dmis_sob2p  fob3  sob3p Dmis_sob3p
	local yvars DDisagree_fob2p DExpDisagree_sob2p Dmis_sob2p    DDisagree_fob3p DExpDisagree_sob3p Dmis_sob3p

* Initial control variables (diagnostics)
local xvars female c1 c2_2_rp c3 c5 f6 age_dif  level_1 level_2 edu_WmoreH  inactiveM unemployedM estrato_12 estrato_3 tec_2 


		foreach sample in `samples' {

********************************************************
* Diagnostics for attrition (who is in sample3_ext)
********************************************************

reghdfe `sample' D `xvars' `bel2', absorb(`strat') vce(cluster id)
	/* sign : D, female, age, level_1, inactivM, tec_2 */
reg sample2 sample1_ext `xvars' `bel2' `strat' if D==0,   vce(cluster id)
	/* sign: being in whatsapp; (W more edu H) ; sob2p; tec_2 */	
reg sample2 sample1_ext `xvars' `bel2' `strat' if D==1,   vce(cluster id)
	/* sign: being in whatsapp;female; age; age_dif;  tec_2; sob2m: inactiveF */

***********	
reg sample2 sample1_ext##D `xvars' `bel2' `strat' ,   vce(cluster id)
	/* sign: being in whatsapp;female; age;  tec_2 */
***********		
	
reg sample1_ext D `xvars' `bel2' `strat' ,   vce(cluster id)
	/* sign: female c5 age, age dif, (level_1), inactiveM, tec_2 */

*********************************************************	
* Diagnostics for selective exposure (who takes D)
*********************************************************

reghdfe D `xvars' `bel2' if `sample' == 1, absorb(`strat') vce(cluster id)

reg D `xvars' `bel2' `strat' if sample3 == 1,   vce(cluster id)
	/* sign: sob2p */
	
****************************************************
* According to diagnostics, our working controls:
* c5 level_1 tec_2 fob2 sob2p
* Working control set for main outcome models
local vars_sample3 `xvars' `bel2'

********************************************************************************
* IPW weights for selection (attrition weights for being in `sample')
********************************************************************************
ta sample3 D ,m
ta sample3 D if sample3_ext ==1,m

cap drop ps_attrition 
cap drop w_*
cap drop wa_*
cap drop ws_*

	probit `sample' D `xvars' `bels' `strat', vce(cluster id)
	*probit `sample' D `xvars' `bels' `strat', vce(cluster id)

predict ps_attrition, pr 

gen w_`sample' = 1 / ps_attrition
quietly summ w_`sample'
local p_uncond = r(mean)
gen ws_`sample' = `p_uncond' / ps_attrition   // normalized (if needed)
gen wa_`sample' = 1                           // unweighted

* Default selection weight (you can switch to ws_ if you prefer normalized)
local weight w_`sample'
summ `weight' 
bys D : su `weight' if sample3_ext ==1
bys D : su `weight' ws_`sample' ps_attrition if sample3_ext ==1
local reps $reps

********************************************************************************
* STORAGE MACROS FOR THE THREE TABLES
********************************************************************************
eststo clear
local j = 1

local estlist_ols_base ""   // 7A: Baseline OLS (no weights) – one model per outcome
local mtitles_ols_base ""

local estlist_ols_ipw ""    // 7B: OLS + selection IPW – one model per outcome
local mtitles_ols_ipw ""

local estlist_ipwra ""      // 7C: IPWRA ATT/PO + OLS-IPW D – one model per outcome
local mtitles_ipwra ""


********************************************************************************
* LOOP OVER OUTCOMES
********************************************************************************

foreach y of local yvars{

    /********************************************************************
     * 1. BASELINE OLS (no selection weights): pooled model only
     *    Male/Female enter as stats rows
     ********************************************************************/

    *** (1.1) Pooled OLS – All (unweighted)
    reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1, absorb(`strat') vce(cluster id)
    scalar N_all = e(N)
    eststo OLS_base_`j'

    *** Means of Y by subgroup (unweighted)
    quietly summarize `y'_w3 if `sample' == 1
    scalar ymean_all  = r(mean)

    quietly summarize `y'_w3 if `sample' == 1 & female == 0
    scalar ymean_male = r(mean)

    quietly summarize `y'_w3 if `sample' == 1 & female == 1
    scalar ymean_fem  = r(mean)

    *** (1.2) OLS – Male only (unweighted)
    reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1 & female == 0, ///
        absorb(`strat') vce(cluster id)
    scalar N_male     = e(N)
    scalar b_D_male   = _b[D]
    scalar se_D_male  = _se[D]
    scalar t_D_male   = b_D_male / se_D_male
    scalar p_D_male   = 2*(1 - normal(abs(t_D_male)))

    *** (1.3) OLS – Female only (unweighted)
    reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1 & female == 1, ///
        absorb(`strat') vce(cluster id)
    scalar N_fem      = e(N)
    scalar b_D_fem    = _b[D]
    scalar se_D_fem   = _se[D]
    scalar t_D_fem    = b_D_fem / se_D_fem
    scalar p_D_fem    = 2*(1 - normal(abs(t_D_fem)))

    *** RITEST p-values for D (All, Male, Female), with capture to avoid r(111)
    * All
    capture noisily ritest D _b[D], ///
        cluster(id) strata(`strat') reps(`reps') : ///
        reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1, ///
        absorb(`strat') vce(cluster id)
    scalar p_ri_all = .
    if (_rc == 0) {
        matrix P_all = r(p)
        scalar p_ri_all = round(P_all[1,1], .001)
    }

    * Male
    capture noisily ritest D _b[D], ///
        cluster(id) strata(`strat') reps(`reps') : ///
        reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1 & female == 0, ///
        absorb(`strat') vce(cluster id)
    scalar p_ri_male = .
    if (_rc == 0) {
        matrix P_m = r(p)
        scalar p_ri_male = round(P_m[1,1], .001)
    }

    * Female
    capture noisily ritest D _b[D], ///
        cluster(id) strata(`strat') reps(`reps') : ///
        reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1 & female == 1, ///
        absorb(`strat') vce(cluster id)
    scalar p_ri_fem = .
    if (_rc == 0) {
        matrix P_f = r(p)
        scalar p_ri_fem = round(P_f[1,1], .001)
    }

    *** Attach all scalars to the pooled OLS model for 6A
    estimates restore OLS_base_`j'
    estadd scalar N_all      = N_all,      replace
    estadd scalar N_male     = N_male,     replace
    estadd scalar N_fem      = N_fem,      replace

    estadd scalar ymean_all  = ymean_all,  replace
    estadd scalar ymean_male = ymean_male, replace
    estadd scalar ymean_fem  = ymean_fem,  replace

    estadd scalar b_D_male   = b_D_male,   replace
    estadd scalar se_D_male  = se_D_male,  replace
    estadd scalar p_D_male   = p_D_male,   replace

    estadd scalar b_D_fem    = b_D_fem,    replace
    estadd scalar se_D_fem   = se_D_fem,   replace
    estadd scalar p_D_fem    = p_D_fem,    replace

    estadd scalar p_ri_all   = p_ri_all,   replace
    estadd scalar p_ri_male  = p_ri_male,  replace
    estadd scalar p_ri_fem   = p_ri_fem,   replace

    eststo OLS_base_`j'

    /********************************************************************
     * 2. OLS WITH SELECTION IPW WEIGHTS: pooled model only
     ********************************************************************/

    *** (2.1) Pooled weighted OLS – All
    reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1 [pweight = `weight'], ///
        absorb(`strat') vce(cluster id)
    scalar N_all_w = e(N)
    eststo OLS_ipw_`j'

    *** Weighted means of Y by subgroup
    quietly summarize `y'_w3 [aweight = `weight'] if `sample' == 1
    scalar ymean_all_w  = r(mean)

    quietly summarize `y'_w3 [aweight = `weight'] if `sample' == 1 & female == 0
    scalar ymean_male_w = r(mean)

    quietly summarize `y'_w3 [aweight = `weight'] if `sample' == 1 & female == 1
    scalar ymean_fem_w  = r(mean)

    *** (2.2) Weighted OLS – Male
    reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1 & female == 0 [pweight = `weight'], ///
        absorb(`strat') vce(cluster id)
    scalar N_male_w     = e(N)
    scalar b_D_male_w   = _b[D]
    scalar se_D_male_w  = _se[D]
    scalar t_D_male_w   = b_D_male_w / se_D_male_w
    scalar p_D_male_w   = 2*(1 - normal(abs(t_D_male_w)))

    *** (2.3) Weighted OLS – Female
    reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1 & female == 1 [pweight = `weight'], ///
        absorb(`strat') vce(cluster id)
    scalar N_fem_w      = e(N)
    scalar b_D_fem_w    = _b[D]
    scalar se_D_fem_w   = _se[D]
    scalar t_D_fem_w    = b_D_fem_w / se_D_fem_w
    scalar p_D_fem_w    = 2*(1 - normal(abs(t_D_fem_w)))

    *** Attach scalars to pooled weighted OLS model for 6B
    estimates restore OLS_ipw_`j'
    estadd scalar N_all      = N_all_w,      replace
    estadd scalar N_male     = N_male_w,     replace
    estadd scalar N_fem      = N_fem_w,      replace

    estadd scalar ymean_all  = ymean_all_w,  replace
    estadd scalar ymean_male = ymean_male_w, replace
    estadd scalar ymean_fem  = ymean_fem_w,  replace

    estadd scalar b_D_male   = b_D_male_w,   replace
    estadd scalar se_D_male  = se_D_male_w,  replace
    estadd scalar p_D_male   = p_D_male_w,   replace

    estadd scalar b_D_fem    = b_D_fem_w,    replace
    estadd scalar se_D_fem   = se_D_fem_w,   replace
    estadd scalar p_D_fem    = p_D_fem_w,    replace

    eststo OLS_ipw_`j'

    /********************************************************************
     * 3. IPWRA – ATT/PO for All, Male, Female (analytic SEs)
     *    Use selection weights as pweights; attach to pooled weighted model
     ********************************************************************/

    *** (3.1) IPWRA – All
    teffects ipwra (`y'_w3 `y' `strat') ///
                   (D `vars_sample3' `strat') ///
                   if `sample' == 1 [pweight = `weight'], ///
                   atet vce(cluster id)

    matrix bA = e(b)
    matrix VA = e(V)

    scalar ATET_all    = bA[1,1]
    scalar ATET_se_all = sqrt(VA[1,1])
    scalar ATET_pv_all = 2*(1 - normal(abs(ATET_all/ATET_se_all)))

    scalar PO_all      = bA[1,2]
    scalar PO_se_all   = sqrt(VA[2,2])
    scalar PO_pv_all   = 2*(1 - normal(abs(PO_all/PO_se_all)))

    *** (3.2) IPWRA – Male
    teffects ipwra (`y'_w3 `y' `strat') ///
                   (D `vars_sample3' `strat') ///
                   if `sample' == 1 & female == 0 [pweight = `weight'], ///
                   atet vce(cluster id)

    matrix bA = e(b)
    matrix VA = e(V)

    scalar ATET_male    = bA[1,1]
    scalar ATET_se_male = sqrt(VA[1,1])
    scalar ATET_pv_male = 2*(1 - normal(abs(ATET_male/ATET_se_male)))

    scalar PO_male      = bA[1,2]
    scalar PO_se_male   = sqrt(VA[2,2])
    scalar PO_pv_male   = 2*(1 - normal(abs(PO_male/PO_se_male)))

    *** (3.3) IPWRA – Female
    teffects ipwra (`y'_w3 `y' `strat') ///
                   (D `vars_sample3' `strat') ///
                   if `sample' == 1 & female == 1 [pweight = `weight'], ///
                   atet vce(cluster id)

    matrix bA = e(b)
    matrix VA = e(V)

    scalar ATET_fem    = bA[1,1]
    scalar ATET_se_fem = sqrt(VA[1,1])
    scalar ATET_pv_fem = 2*(1 - normal(abs(ATET_fem/ATET_se_fem)))

    scalar PO_fem      = bA[1,2]
    scalar PO_se_fem   = sqrt(VA[2,2])
    scalar PO_pv_fem   = 2*(1 - normal(abs(PO_fem/PO_se_fem)))

    *** Create IPWRA container: copy pooled weighted OLS and attach ATT/PO scalars
    estimates restore OLS_ipw_`j'
    eststo IPWRA_`j'

    estimates restore IPWRA_`j'
    estadd scalar ATET_all    = ATET_all,    replace
    estadd scalar ATET_se_all = ATET_se_all, replace
    estadd scalar ATET_pv_all = ATET_pv_all, replace

    estadd scalar PO_all      = PO_all,      replace
    estadd scalar PO_se_all   = PO_se_all,   replace
    estadd scalar PO_pv_all   = PO_pv_all,   replace

    estadd scalar ATET_male     = ATET_male,     replace
    estadd scalar ATET_se_male  = ATET_se_male,  replace
    estadd scalar ATET_pv_male  = ATET_pv_male,  replace

    estadd scalar PO_male       = PO_male,       replace
    estadd scalar PO_se_male    = PO_se_male,    replace
    estadd scalar PO_pv_male    = PO_pv_male,    replace

    estadd scalar ATET_fem      = ATET_fem,      replace
    estadd scalar ATET_se_fem   = ATET_se_fem,   replace
    estadd scalar ATET_pv_fem   = ATET_pv_fem,   replace

    estadd scalar PO_fem        = PO_fem,        replace
    estadd scalar PO_se_fem     = PO_se_fem,     replace
    estadd scalar PO_pv_fem     = PO_pv_fem,     replace

    eststo IPWRA_`j'

    ********************************************************************
    * Append models and column titles
    ********************************************************************
    local ylab : variable label `y'_w3
    if "`ylab'" == "" local ylab "`y'_w3"

    local estlist_ols_base "`estlist_ols_base' OLS_base_`j'"
    local mtitles_ols_base `"`mtitles_ols_base' "`ylab'""'

    local estlist_ols_ipw  "`estlist_ols_ipw'  OLS_ipw_`j'"
    local mtitles_ols_ipw  `"`mtitles_ols_ipw'  "`ylab'""'

    local estlist_ipwra    "`estlist_ipwra'    IPWRA_`j'"
    local mtitles_ipwra    `"`mtitles_ipwra'    "`ylab'""'

    local j = `j' + 1
}


*=====================================================
* Table 7A – BASELINE OLS (no selection weights)
*=====================================================
#delimit ;
local opts_ols_base style(tex)  
  varlabels(_cons Constant D "Treatment" female "Female (=1)") 
  cells(b(star fmt(%9.3f)) se(par)) 
  stats(N_all N_male N_fem 
        ymean_all ymean_male ymean_fem
        b_D_male se_D_male p_D_male
        b_D_fem  se_D_fem  p_D_fem
        p_ri_all p_ri_male p_ri_fem,
        fmt(%9.3f)
        label("N (All)" 
              "N (Male)"
              "N (Female)"
              "Mean Dep. Var (All)"
              "Mean Dep. Var (Male)"
              "Mean Dep. Var (Female)"
              "OLS β(D) Male"    "OLS SE Male"    "OLS p-val Male"
              "OLS β(D) Female"  "OLS SE Female"  "OLS p-val Female"
              "RI p-val D (All)" "RI p-val D (Male)" "RI p-val D (Female)")) 
  label collabels(none) msign(--) nolz varwidth(1) modelwidth(14) 
  starlevels(* 0.10 ** 0.05 *** 0.01)
;
#delimit cr

estout `estlist_ols_base' using "$outdir/Table7_SpousalBeliefs_y2_OLS_baseline_`sample'.tex", ///
    replace `opts_ols_base' mlabels(`mtitles_ols_base')

*=====================================================
* Table 7B – OLS WITH SELECTION IPW WEIGHTS
*=====================================================
#delimit ;
local opts_ols_ipw style(tex)  
  varlabels(_cons Constant D "Treatment" female "Female (=1)") 
  cells(b(star fmt(%9.3f)) se(par)) 
  stats(N_all N_male N_fem 
        ymean_all ymean_male ymean_fem
        b_D_male se_D_male p_D_male
        b_D_fem  se_D_fem  p_D_fem,
        fmt(%9.3f)
        label("N (All)" 
              "N (Male)"
              "N (Female)"
              "Mean Dep. Var (All, weighted)"
              "Mean Dep. Var (Male, weighted)"
              "Mean Dep. Var (Female, weighted)"
              "OLS-IPW β(D) Male"    "OLS-IPW SE Male"    "OLS-IPW p-val Male"
              "OLS-IPW β(D) Female"  "OLS-IPW SE Female"  "OLS-IPW p-val Female")) 
  label collabels(none) msign(--) nolz varwidth(1) modelwidth(14) 
  starlevels(* 0.10 ** 0.05 *** 0.01)
;
#delimit cr

estout `estlist_ols_ipw' using "$outdir/Table7_SpousalBeliefs_y2_OLS_IPW_`sample'.tex", ///
    replace `opts_ols_ipw' mlabels(`mtitles_ols_ipw')

*=====================================================
* Table 7C – IPWRA (analytic SEs), only D in coefficient panel
*=====================================================
#delimit ;
local opts_ipwra style(tex)  
  cells(b(star fmt(%9.3f)) se(par)) 
  stats(ATET_all ATET_se_all ATET_pv_all 
        PO_all   PO_se_all   PO_pv_all
        ATET_male ATET_se_male ATET_pv_male
        PO_male   PO_se_male   PO_pv_male
        ATET_fem  ATET_se_fem  ATET_pv_fem
        PO_fem    PO_se_fem    PO_pv_fem,
        fmt(%9.3f)
        label("ATT All"      "Std. Err. ATT All"   "P-value ATT All"
              "PO All"       "Std. Err. PO All"    "P-value PO All"
              "ATT Male"     "Std. Err. ATT Male"  "P-value ATT Male"
              "PO Male"      "Std. Err. PO Male"   "P-value PO Male"
              "ATT Female"   "Std. Err. ATT Female""P-value ATT Female"
              "PO Female"    "Std. Err. PO Female" "P-value PO Female")) 
  label collabels(none) msign(--) nolz varwidth(1) modelwidth(14) 
  starlevels(* 0.10 ** 0.05 *** 0.01)
;
#delimit cr

estout `estlist_ipwra' using "$outdir/Table7_SpousalBeliefs_y2_IPWRA_main_`sample'.tex", ///
    replace `opts_ipwra' mlabels(`mtitles_ipwra') ///
    keep(D) varlabels(D "Treatment") 

}